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In this paper a feasible method of implementing the 
Maximum Principle is given based on an adaptive random 
search algorithm which utilizes direct measurements of 
the boundary cost -function hyper -surfaces . No restric- 
tions are placed on the continuity of the surfaces or 
the number of valleys. The adaption enhances convergence 
by varying the mean and variance of a probability distribu- 
tion as a function of the past performance. The algorithm 
has both local and global search properties so that "hanging 
up" in local valleys is avoided. 

A hybrid computer implementation of the algorithm is 
discussed. The usefulness of the algorithm is investigated 
experimentally for a fifth-order nonlinear minimum-fuel 
orbit -transfer problem. Convergence is generally obtained 
within several thousand iterations (one to two minutes) . 
Details are given on the manner in which the system iterates. 


typical solutions obtained for a wide range of situations, 
and the convergence properties of several variations of 
the basic algorithm. Cross sections through the boundary 
hyper -surfaces reveal the striking irregularities which 
high-order systems can have and, hence, demonstrate the 
effectiveness of the adaptive random search approach in 
coping with them. 

INTRODUCTION 

The Pontryagin Maximum Principle is an exceedingly 
elegant and powerful theory for determining the optimal 
control of dynamic systems describable by nonlinear dif- 
ferential equations with bounded control. Although there 
has been a flurry of activity for several years in its 
application to low-order systems, there have been few 
applications to high -order systems. For high-order systems 
the Maximum Principle yields much information about the 
nature of the solution, but the actual solution is generally 
difficult to obtain. The reason for this is that a diffi- 
cult mixed boundary -value problem is invariably encountered, 
one in which the known boundary conditions are divided between 
initial and terminal values. Furthermore, the mapping from 
the initial to terminal boundary values is not generally 
explicitly known. 
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One possible approach to the boundary -value problem is 
based on some form of the gradient method. However, there 
are several disadvantages in this approach: only local prop- 

erties are utilized, only one local minimum is implicitly 
assumed, and the gradient is not analytically available. 

These difficulties were shown in reference 1 to lead to con- 
vergence problems. Thus the approach is of limited value when 
the surface being studied is multi-peaked, discontinuous, or 
very flat in certain regions. Random search methods would 
seem to be promising in overcoming some of the above deficien- 
cies. Such methods have been used mostly in connection with 
direct search problems (refs. 2 and 3 ). However, in a recent 
study (ref. 4) a fixed-step, random-sign-type search was used 
to implement the Maximum Principle, and application was made 
to linear second- and third-order systems. 

In this paper we will investigate a more general approach 
based on an adaptive random search method for solving the two- 
point boundary -value problem involved in the Maximum Principle. 
The search has both local and global properties. The method 
does not depend on the continuity of the surface being searched 
or the number of minima, and does not require that an analytical 
expression for the surface be known. In the sections to 
follow, the adaptive random search will be described, the 
hybrid computer realization discussed, and the results of 
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an application of the method to a fifth-order nonlinear 
orbit -transfer problem presented. 


THE ADAPTIVE RANDOM SEARCH ALGORITHM 
General Approach 

The class of problems considered here are restricted to 
those for which the Maximum Principle is applicable. Famil- 
iarity with the Maximum Principle is assumed (refs. 5 , 6 , 7 ). 
Thus, we will not review the derivation, the theorems, the 
assumptions involved, or the boundary -value theory. For 
purposes of notation, however, a few remarks are appropriate. 
The system to be controlled is defined by the vector equation. 

x = f ( x , u, t) (l) 

where x = (x^., x 2 , . • • x n ) , u = (uj., u 2 , • • • u m ) , and 
ueU, the control region. Interest will center on fixed-time 
problems because of- its convenience in computer operation. 

It will be desired to take the system from a given state 
x(0) to a final target set S so as to minimize the 


generalized cost function 
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where x Q (t) is the auxiliary state associated with the 
quantity to be minimized (ref. 7) • Some of the possible 
target sets of interest are (a) SeR n , (b) S <= R n , and 
(c) S = X|.eR n , corresponding to a free point, a subset of 
the whole space R n , and a fixed point. The solution to 
the problem invariably requires the solution to the set of 
equations 


u = u(x, p, t) 

x = f(x, u, t) 

P * g(p* X, u, t) 


(3) 


where the known boundary values are divided between the 
initial and terminal values. The final boundary condition 
required on p is dependent on the desired final boundary 
condition on the states. 

The general approach to be used here to obtain 
explicit solutions to the Maximum Principle is based on an 
adaptive random search procedure. A hybrid computer diagram 
of the approach is illustrated in figure 1. We will be 
concerned with the fixed time interval 0 to T; however, 
it is trivial to vary this time in the computer implementa- 
tion. The left half of the figure, indicated by the analog 
computation, is an implementation of the set of equations (3) • 
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The right half, indicated by the digital computation, is the 
algorithm which will enable the boundary conditions involved 
in the Maximum Principle to be satisfied. 

Let us discuss the algorithm in detail. The basic notion 
in the random search algorithm is to select the initial con- 
dition vector for the adjoint equations from a noise source, 
in this case a gaussian distribution with fixed mean m and 
fixed standard deviation a. Conceptually, the mean and stan- 
dard deviation can be separated as indicated in the figure by 
adding the mean m to the output of a gaussian noise source 
with standard deviation a and zero mean. Thus, on any kth 
iteration the vector p is 

p k = m k + | k (4) 

The gaussian source generates a purely random, gaussian, 

n-dimensional vector sequence with zero mean and 

independent components: 


m s E(| k ) = 0 



cov(|j, | k ) = Q 

if J 

= k 

0 

if j 

i k 


where i - (ii> • • • ? n ) 

Q = a 2 I 

and I is the identity matrix. 

Continuing around the loop in figure 1, the value of p k 
as determined by equation (4) becomes p(0), the initial 
condition for the adjoint equation. Since the value of p(0) 
together with the given x(0) is sufficient to define a 



solution of the set (3)> these equations, as represented by 
the left half of figure 1, can be integrated to the terminal 
time T. The final state x^(T) actually achieved will 
generally fail to satisfy the desired terminal state as 
defined by the target set. Similarly, the final values of 
the adjoint variables p^(T) will fail to satisfy the bound- 
ary conditions required depending on the target set. For 
this reason we will introduce a function to measure the dif- 
ference between the actual boundary values and the desired 
boundary values. The deficiencies pointed out in reference 8 
of the scalar -valued metric will be avoided. We will par- 
tially order systems by a vector -valued metric. Since there 
will generally be three distinct types of boundary conditions 
to satisfy at the terminal time, we will take the vector 
metric of the form 

J = Jyw Jp) (6) 

where the components refer to displacement, velocity, and 
adjoint variables, respectively . For two systems S and S' 
we will say S > S' if and only if J > J' , that is, 

Jp > Jp' , Jy > Jy* , and Jp > Jp' . With this concept one is 
concerned with choosing from a certain set of systems a 
noninferior system rather than an optimum system. This 
corresponds more closely with a realistic objective as pointed 
out in reference 8. 
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Assume now that the search is purely random with no 
adaptive characteristics. In this case the dashed lines 
indicated in figure 1 would be absent. Then as a result of 
the random vector sequence the vector sequence 

is generated. Since the search is purely random, the 
algorithm logic decides at each iteration whether the value 
of J has been reduced to zero. 

Satisfying the boundary conditions is slightly more 
involved than the above because in an experimental study it 
is not likely the boundary conditions for x and p will 
ever be precisely achieved. For this reason we will enlarge 
the target set by some small amount and allow the system to 
terminate at any point in the enlarged set. This view is 
more realistic, since there is no reason to demand that a 
practical system satisfy the boundary condition exactly. 

To accomplish this in the random search method, we will 
require of the vector metric 

J < e (7) 

where 

e = ( e D> e v> € p) 

The values of £p and ey are dictated by how closely it is 
desired that the system states approach the desired target 
set. The value of ep is more difficult to choose because 
of its lack of physical interpretation. Fortunately., the 
hybrid computer approach makes it easy to choose a suffi- 
ciently small value by experimentally observing the sensitivity 
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of the solutions to small changes in ep. In an experimental 
study to be (described, later, the solutions were quite insen- 
sitive to variations in ep over a wide range. 


The Adaptive Algorithm 

The pure random search described in the preceding section 
is generally unsatisfactory because of the excessively long 
convergence times, as will be seen in a later example. We 
attempt to improve the convergence properties by making the 
system adaptive by varying the mean m and standard devia- 
tion a as a function of the system's past performance. 

Since performance is determined by the input and output 
sequences, {| k } and {j k }, we will make the mean and standard 
deviation adaptive of the form 





CT k+i =, f^d 1 , i 2 , . . . £ k , J°, J 1 , . . . J k ) 

Some reasonable and easily implemented forms will be discussed 
in the following paragraphs. A fundamental notion in the 
algorithm will be that of a "success" or "failure" defined by 
a function of the cost J on any kth iteration and the 
smallest preceding cost obtained on the 2th iteration. 

The most frequently used definitions were simply: 

success: - J 11 > 0 

i k (9) 

failure: j‘ - J 5 < 0 
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Obviously a successful iteration is both necessary and suffi- 
cient for the system to be noninferior. It might be noted 
that v T hen any component of is less than its corresponding 

£ component , we will not require further reduction for a 
success as long as it remains less than e. 

Implementation of the adaptive part of the algorithm is 
illustrated in figure 1 by the dotted lines and the algorithm 
logic block. The algorithm logic block performs the computa- 
tions in the above equations by utilizing: (l) the p^ 

information from the memory Mp, and (2) the system oerformance 
information, J , from the memory Mj . The output of the 
algorithm logic block is then the mean and standard deviation 
of the distribution for the next iteration as indicated. 

The rationale behind the particular adaptive law to be 
used here is based on the notion of a creeping, expanding, 
and contracting search. The creeping character of the search 
is provided by varying the mean of the distribution so as to 
equal the initial condition of the adjoint vector on the last 
successful iteration. The expansion and contraction charac- 
ter of the search is provided by varying the standard deviation 
such that the search is localized when successful but gradually 
expanded when not successful. These characteristics give the 
algorithm some useful search properties. 
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The adaptive law for the mean value of the distribution 
was taken to be: 


m k+1 = p k 

if 

J k < J^ 1 1 

(10) 

m k 

•if 

j k > j k_i ! 


where the initial values are 

mi 

= 0 and J° 

is the value 


based on the initial values x(0) . That this law is of the 
form (8) can be seen by combining (4) and (10) sequentially. 
Thus, the first few terms of the sequence are 


mi = 0 


m 2 = 


5 1 

0 


if Ji < J° 
if J 1 > J° 


= f/Cft 1 , J°, J 1 ) 


m 3 =( 


|i + | 2 if J 2 < J 1 < J° 

I 1 if J 2 > J 1 < J° 

if J 1 > J°, J 2 < J° 




) « f/U 1 , t 2 , J°, J 1 , J 2 ) 


r 

10 


if J 2 , J 1 > J' 


(ID 


The adaptive law for the standard deviation of the dis- 
tribution, cr, was implemented so as to expand the size of the 
search depending on performance. Let be the metric 

obtained on the last successful iteration. Then the law for 
the standard deviation on any kth iteration is 



I 
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- Cl 

if 

k - 

q < l 



= ct 2 

if 

•k > 

. . . A’ > J 1 

1 

1 

(12) 

• 0*y 

if 

A 

. . . JMr-ih > J 1 J 




where the o^'s are constants such that ax < a 2 < . . . a ^ . 

In other words, is used if there was a success in the last 
q iterations, a 2 if there was no success in the last q iter- 
ations, as if no success in the last 2q iterations, etc. 

These equations are, of course, a simple form of equations (8) . 
Occasionally, depending on the difficulty of the particular 
problem being studied, no success will be achieved in the 7 q 
iterations. For this reason, we repeat the search in equa- 
tions (12) a number of times C (usually twice), and then 
reinitialize the entire search if no success is achieved. 

The parameters which are free to choose are q, 7, and the 
cr^'s; in the example to be discussed later values of q =• 100 
and 7 = 10 were chosen and they do not seem to be critical. 

The choice of a values is more important because it affects 
both the smallness and the largeness of the search. However, 
it is not difficult to choose reasonable values. A reasonable 
lowest value, oi, can be determined by observing on the dis- 
play system (yet to be described) the metric J on every 
iteration, and choosing a value small enough so that the J 
generally varies only slightly. For the orbit -transfer 
problem discussed later, a\ =0.1 volt proved satisfactory. 
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The upper value is chosen to cover some sizable portion of 
the entire space; a value of 10 volts seemed reasonable for 
the 100 volt space available on the analog computer. In 
between, the steps are a geometric progression to enable 
the search to expand rapidly. 

Several additional adaptive strategies were incor- 
porated into the algorithm to provide better convergence 
properties and greater versatility in certain situations. 

They are itemized by the following: 

1. Single-step strategy -- This strategy is intended 
to take advantage of favorable local properties of 
the surface. It is based on the notion that the 
step following a successful step k should not be 
random but deterministic in the same direction and 

k+i 

of the same amount. That is, p is defined by 

pk+i _ m k+i + |k if jk-i _ jk > o (13) 

The cost in time, one iteration, is insignificant, 
and the benefits are substantial as will be seen in 
a later example. 

2. Threshold strategy -- In this strategy the require- 

ment for a success is modified so that the differ- 
ence between the cost for the kth iteration and 
that for the last successful iteration (i.e., the 
left side of equations ( 9 )) must exceed some thresh- 
old value dependent on the cost function. We take 
the case in which the threshold is simply where 

0 < q <1. Thus, a success is defined by 
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J l - J k > ri J 1 (14) 

This strategy has the desirable characteristic that 
smaller improvements are required as the minimum is 
approached. 

3 . Localized end -search strategy -- This end-search 

strategy is intended to localize the search when in 
the immediate neighborhood of the minimum. This is 
done by reducing the standard deviation in equa- 
tions (l2) so that a varies from clg± to a cr^ 
k 

when J <3, where 0 < a < 1, and p is some small 
value generally on the order of 2e . This strategy 
was used sparingly since it was beneficial in reducing 
convergence time only under certain conditions. 
Although a gradient method could be used in this final 
phase y the random search method accomplishes the same 
objective without an implementation change. 

k* Initial search --It was found advantageous to control 
the distribution from which the first successful 
adjoint vector p(0) is obtained. This was done by 
starting the search with a square distribution of 
variable width, -W to +W. When a success is obtained 
the search continues with a gaussian distribution as 
described above. The rationale behind the square 
distribution is that since no information is available 


on the best starting value, the weighting within the 
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limits should be equal. The size of W is a com- 
promise: It must be large enough to include possible 

solutions but not so large that the volume being searched 
is excessive. The choice of W will be dependent, of 
course, on the particular problem; it is not difficult 
to choose experimentally a reasonable value. Values of 
15 to 30 volts are typical for the later example 
problem. 

It is clear the behavior of the adaptive algorithm is 
different from the pure random or nonadaptive search described 
before. Now the vector metric J which measures the disparity 
between desired and actual terminal boundary conditions is 
sequentially reduced rather than being reduced in one iteration. 
The way in which this occurs is illustrated in figure 2, where 
a representative surface is given in only two dimensions. Due 
to its adaptive character, after any success the mean of the 
distribution moves to the last successful p. At this point 
the search starts locally and increases in a geometric pro- 
gression until the next success is obtained. In this way, 
the successive values of J may jump from one valley to 
another as indicated by the numbered points until a J less 
than the required e is reached. The algorithm logic block 
decides when this condition occurs. 

The virtues of the random search approach are clear. It 
has desirable local and global search properties so that "hanging 
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up t! in local valleys as with gradient methods is avoided. 
Further, it will not matter if the surface is discontinuous or 
how many peaks or valleys there are. The main question is, 
of course, the convergence time. This is best studied 
experimentally on an example to be discussed in a later 
section. 

IMPLEMENTAT ION 

Considerations in the Choice of Computer System 

The hybrid computer proved to be the most feasible way 
in which to implement the random search algorithm. The reasons 
for this and a comparison with the alternatives are worth 
discussing. 

In general there are three basic computational techniques 
available to the experimenter for the implementation: analog, 
digital, or a combination of the two, hybrid. Of prime 
importance in the selection of the computer system is a con- 
sideration of the number of iterations required to find a 
solution, and this number is not known a priori. For a 
second-order system, pilot studies indicated something on 
the order of 10 2 iterations to obtain a solution. Since 
the volume of the space increases so rapidly with the order 
of the system, one might expect several thousand iterations 
to be necessary for an increase in system order to perhaps 
five or six. Thus, we see that the time per iteration 
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will be of critical importance and will largely dictate the 
means for implementing the search algorithm. 

Each iteration can be divided into two steps: 

(l) integration of the equations of motion on the interval 
[0, T], and (2) execution of the algorithm. For the orbit - 
transfer problem to be discussed later, an IBM 709*+ machine 
requires 1-10 seconds to perform the integration in step 1. 
An analog computer, however, performs the same integration 
in 1-10 milliseconds with an accuracy of about 5 percent 
relative to the digital machine solutions. Thus for this 

o 

specific example, the analog computer is about 10 faster 
than the digital computer in performing the integration 
in step 1. The second step is best accomplished digitally. 
The time to perform the second step on a digital computer 
of speed comparable to the IBM 709*+ is the same order of 
magnitude as the time required for the analog to perform 
step 1. Therefore, a great saving in computer time can be 
realized over a completely digital simulation by a hybrid 
approach. It is worth noting that an alternative approach 
was investigated utilizing pseudo-hybrid techniques, that 
is, an analog computer and something less than a digital 
computer. However, our experience shows that inaccuracies, 
limited storage, and limited flexibilities in logical 
operations seriously limit the feasibility of this approach. 
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The Hybrid System 

In the hybrid implementation the analog computer was 
delegated the task of solving the state, adjoint, and control 
equations, as given in equations (3). It also served as the 
point at which the operator exercised manual control over 
the hybrid system. The digital computer was required to 
calculate the metric, provide storage, implement the algorithm, 
generate the initial conditions for the adjoint equations, 
and finally, oversee the sequencing of events of the iterate 
cycle. 

Figure 3 is a hardware diagram of the hybrid system used. 
Shown are the two basic elements of the simulation, the 
analog and digital computers along with their coupling system, 
and peripherals. The coupling system is comprised of two 
distinct parts: (a) the Linkage System and (b) the Control 

Interface System. A discussion of the hardware used in these 
subsystems is given in the four sections to follow. A final 
section discusses the sequencing of events through the sub- 
systems during one iteration cycle in order to better describe 
the functioning of the hybrid system as a whole. 

A. Digital Computer — The digital computer used in 
the optimization program was an Electronic Associates, Inc. 
(EAI) model 8400. It is a medium-sized, high-speed computer 
which is designed to operate in a hybrid atmosphere. The 
particular machine used has 16,000 words of core memory 
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with 32 bits per word. Memory cycle time is two "microseconds 
The machine uses parallel operation for maximum speed. Float 
ing point operations are hardware implemented. Programming 
languages available include a MACRO ASSEMBLY language and 
FORTRAN IV. The optimization program was coded in MACRO 
ASSEMBLY in order to keep the execution time to a minimum. 

The instruction repertoire includes special commands by which 
discrete signals can be sent to or received from the external 
world. External interrupts are provided which can trap the 
computer to a specific cell in memory. In an example to be 
discussed later the optimization program utilized about 
8,000 words of storage. Of this about 1,000 comprise the 
actual optimization executive program, the remaining "J ,000 
being used for subroutines, monitor and on-line debugging 
and program modification routines. 

The peripherals of the digital computer include mag- 
netic tapes, card reader, and printer. The two former 
devices were used for program input and storage while the 
latter device was used for data logging. 

B. Analog Computer — The analog hardware consisted 
of an Electronic Associates 231R-V analog computer. The 
state equations, adjoint equations, and the control logic 
were programmed in standard fashion. Consequently, analog 
schematics were thought not essential in this paper. 
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The analog computer serves as the point at which mode 
control of the hybrid computer is accomplished. By manual 
selection of switches either of two modes can be commanded: 
(l) In the "search" mode the analog computer operates in a 
high-speed repetitive manner. Such operation is accomplished 
by controlling the mode of the individual integrators with 
an appropriate discrete signal. This signal is a two -level 
signal which is generated on the control interface in con- 
junction with the digital computer and, depending on the 
level, holds an integrator in either "operate" or "initial 
condition" mode. (2) In the "reset" mode, the integrators 
are placed in their initial -condition mode and held there. 

Particular equipment worth pointing out are the track- 
store units, D/A switches, and comparators with which the 
control logic was implemented. At the end of the operate 
period T, the digital computer reads a number of variables 
essential to its functioning. Because of the high repeti- 
tive speeds used, the value of a variable could change 
considerably between the end time T and a later time when 
it is actually converted. Thus, track-store units were 
used to hold the variables at their respective values at 
time T until the digital read all of the values. The 
control logic requires on-off type switching and the high- 
speed electronic comparators and electronic switches were 
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quit e necessary for proper operation. Nonlinear operations 
such as multiply and divide provided no particular difficul- 
ties , in their normal operation, although a square root 
operation did require the use of a diode function generator. 

For continuous type output, a display console was con- 
nected to the analog computer to provide visual readout of 
variables. The display contained a cathode ray tube (CRT) 
which could simultaneously display up to four channels, 
and enabled photographic records to be taken of the display 
quantities. The display was extremely helpful in determining 
if the algorithm was functioning properly. Other continuous 
analog data useful for determining proper functioning were the 
p^ and ; these could be recorded on a pen recorder since 
their rate of change was low. 

C. Control Interface -- The control interface between 
the analog and digital system is an Electronic Associates, 

Inc. DOS 350 (see fig. 3). It is through this unit that 
the iteration process is controlled. An important task 
allocated to this subsystem is the operate-time control. This 
function is implemented through the use of a counter and 
is the key element in the control of all timing in the hybrid 
simulation. The counter is driven from a high-frequency 
source in the interface system allowing for a very high degree 
of resolution in the simulated operate-time. The counter is 
constructed by patching modular blocks which can be combined 
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to give a wide range of simulated operate-time. Specific 
times within this range are selected with thumb-wheel switches. 
Also, the interface allows the digital computer to use any 
conditions in the analog computer which can be represented 
by discrete variables (binary levels) and to send discrete 
signals to the analog system to be used as control levels or 
indicators. An example of the former would be the hybrid 
system mode control which merely amounted to the operator 
depressing the "reset" or "search" switch on the analog 
console. This action sets a binary level which is then 
sensed by the digital computer. An example of the latter 
situation is when the digital sends the operate command to 
the operate-time counter. The interface system allows 
patching of Boolean functions. Hence, some of the logic 
operations required for timing pulses, event signals, and 
other like operations were very effectively programmed on 
it. 

D. Linkage System -- The linkage system which is 
shown in figure 3 is that part of the system which houses 
the conversion equipment, the A/D and D/A converters. It 
is through here that all of the data passes between the 
analog and digital portions of the simulation. The linkage 
system is controlled by command from the digital computer. 

Input to the digital computer is through the A/D 
converter which has preceding it a channel selection device 



or multiplexer to select the analog channel which is to be 
converted. Conversions were done sequentially through the 
analog channels at a maximum rate of 80,000 samples per 
second from channel to channel. 

Output to the analog used the D/A converters with each 
data channel having its own conversion unit. The maximum 
conversion rate of the D/A’s used is 250,000 conversions 
per second. 

E. Sequencing of Events During One Iteration — The 
sequencing of events during one iteration cycle are depicted 
in figure 4 . The instants of time ti, t 2 , • • • , ts shown 
in this figure are considered fixed relative to each other, 
and ti is conveniently regarded to be the start of the 
iteration cycle. We will consider the cycle to begin at 
ti with the analog integrators in an operate mode. As 
discussed previously, the elapsed time (t 2 - ti) is controlled 
by a counter on the interface system. At t 2 an interrupt 
pulse is generated on the control interface which is sent to 
the digital computer signaling it to commence its operations. 
Simultaneously, the pulse is sent to the analog to instruct 
the track -store units to hold their respective values which 
they possessed at time t 2 * During the interval (t 3 - t 2 ) 
the digital computer reads these analog variables with the 
A/D converter. At t 3 the digital sends a pulse via the 
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interface to the analog console which commands the integrators 
to an initial condition mode. At when the data required 

by the algorithm has been generated the D/A converters send 
these values to the appropriate points in the analog portion 
of the simulation. The digital machine allows enough time 
for the transients to settle in the initial condition circuits 
of the analog before sending a command at ts that places the 
integrators in an operate mode and starts the counter. Since 
ts and ti are the same event, we merely repeat the above 
sequence for repetitive operation. 

Some specific numerical values might be of interest. 

The total iterate time (ts - ti) is primarily composed of 
two parts: (l) (t 2 - t x ) which in a later example problem 

was scaled in the simulation to 2.5 milliseconds, and (2) 

(ts - t 2 ) , which was primarily determined by the speed of 
the digital machine in computing, converting, and generating 
random numbers; this latter period was on the order of 
7-5 milliseconds. Thus, the total iterate time for the above 
situation is on the order of 10 milliseconds (or 100 iterations 
per second) . This figure is dependent on the control problem 
chosen and the exact form of the algorithm implemented. 

The Algorithm Flow Graph 

A program flow graph which shows the operation of the 
iteration process is displayed in figure 5- This basically 



is the flow graph of the algorithm and the iteration control 
sequences utilized by the hybrid system. Note the inclusion of 
the event times ti, . . . discussed earlier in connection 
with figure 4. The program is continuously recycling in a 
high-speed repetitive fashion. 

There are three basic loops which correspond to the 
three system modes in the optimization program: a reset 

loop, a search loop, and an end -state loop. The reset loop 
is for the purpose of initializing the program. The search 
loop is that portion of the program which uses the algorithm 
to search for a solution to the problem. The end-state loop 
is entered by the digital program when a solution is found, 
and is used for the generation of graphic displays. The 
operator manually selects the search or reset mode as dis- 
cussed in the section dealing with the analog computer. 

A. Reset Loop -- The reset loop is the portion of 
the repetitive operation cycle which initializes the program 
and prepares it for the search mode. When the system is in 
reset mode, the integrators of the analog computer remain 
in their initial condition state. The flow of the reset 
cycle is shown in figure 6 where that particular loop is 
emphasized by line weight, figure 6 being the same as 
figure 5 otherwise. It is while the system is in the reset 
mode that the digital program is first entered, the point of 
entrance being designated by START in the figure. The first 
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operation performed is that of initialization. It is during 
this process that all the program variables are set to their 
initial states. Also, the line printer is initialized and 
and the data header printed. Once complete, the interrupt 
line is enabled and the analog computer signaled that an 
iteration cycle is to begin. It is at this point that the 
counter which controls the operate time of the analog is 
started. The digital computer than halts and waits for an 
interrupt to signal that the period T has ended. 

When the interrupt ocfcurs, the digital program proceeds 
to the next operation where the values of the states x^(T) 
and p k (T) of the analog computer are read by the A/D converter. 
Once the inputs to the digital computer have been read, a 
pulse is sent to the analog computer to place the integrators 
in their initial condition mode. Since the integrators are 
already in initial condition state due to the reset mode, 
this pulse has no effect but it is needed later when the 
system is in the search mode or the end-state mode where the 
integrators are not so constrained. 

The next operation performed is that of calculating the 
metric, J^, based on the samples of the state obtained in 
the last block. This value is peculiar to the reset mode 
and, therefore, is designated J° and saved as such. Once 
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the metric is calculated, the digital computer tests the 
system mode and branches to that portion of the program 
which is exclusive to the reset loop. 

The program now initializes the algorithm by setting 
/ : p^(0) =0 and = J°. The values of the components 
of p k (0) are then generated using a uniformly distributed 
noise source, the space of the noise being [ -W, W]. Note 
that in reset mode p k (0) = 

At this point the reset loop returns to the mainstream 
of the program by entering the output portion of a cycle. 

It is here that all of the quantities required at the analog 
computer are converted to analog form by the D/A converters. 
The data sent to the analog include p k (0), p^(0), J^, and 
J k . Note that p k (0) is required by the analog computer 
whereas the others are displayed to determine if the algorithm 
is functioning properly. 

Finally the digital program enables the interrupt line 
and signals the analog computer and control interface to 
begin another cycle. As before, the program now goes into a 
halted state waiting for the interval T to pass. Cycling 
in the reset loop continues until the operator is ready to 
begin a search cycle and does so by putting the hybrid system 
in the search mode. 

B. The Search Loop — It is in the search mode that the 
algorithm is used to seek an optimal solution to the control 



-28- 


problem implemented. Search mode is selected by the operator 
at any time and is considered to begin on the first iteration 
after his doing so. 

During the first iteration of search mode, the analog 
computer solves the system, adjoint and control equations 
using the p(0) T s established in the last cycle through the 
reset loop. The digital computer then receives the time T 
interrupt signal and proceeds to read the states x^(T) and 
P k (T) as shown in figure 7 * Once read, the integrators are 
set back to initial condition where they will await the next 
p(0) to be generated by the algorithm. The next operation 
computes the metric as was done in the reset loop; only this 
time the J k is based on the actual solution of the system 
equations at time T. 

The system mode test is next and results in a branch to 
that portion of the program that solves the algorithm and 
generates the p(0)'s based on the adaptation principles. 

The search branch begins with a. test to determine if the 
system is in end state. Since this is the first time through 
the search loop, end-state condition cannot have been established 
so the program proceeds to the operation which tests the metric. 
Here all of the components of the metric are tested according 
to equations (9) and the condition of a success or failure 
ascertained as discussed before. 
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At this point in the discussion we shall take the results 
of the metric test to be a success so that that portion of 
the search loop can be examined. The first operation in this 
segment up-dates the memory of the algorithm by setting the 
last good values of the p(0)'s to p^(o), the vector which 
gave this success, and by setting the last good metric 
to jk. Next a test is made on the metric to see if it 
meets the requirement of a solution, i.e., J^< e. Assume 
for the moment that it does not and that the program proceeds 
out the lower branch of the test. The program then generates 
the next try for the adjoint initial conditions on the basis 
of the single-step strategy. Following the establishment of 
p k+1 (0) the program tests J k to determine if the end-search 
strategy should be employed by comparing with 5. If 

is less, the noise standard deviation sequence is reduced 
through multiplication by the constant a; if is not less 

than 5, the standard deviation is unaltered. Then after 
setting a to the starting value Oi, the program execution 
returns to that part which is common to all loops, the output 
portion. The program then starts the next iteration exactly 
as it does in the reset loop. 

If the test of the metric had resulted in a failure, the 
program would have branched to the right after the metric test 
shown in the flow graph. Here the primary function is to 
generate a new set of adjoint initial conditions using the 
noise generator as discussed in the description of the 
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adaptive algorithm. Before the actual generation of p k+1 (0) 
can take place, however, certain of the algorithm parameters 
must be examined. These parameters control the maximum time 
at a given a level, the changing of cr levels, and whether 
or not to start the search anew. First, a test is performed 
to find out if q consecutive trials have resulted in 
failures. If not, a remains fixed and the program executes 
the jump ahead shown; if true, a is incremented to the next 
value in the sequence (<Ji, . . ., Oy) . Next is a test to 
determine if a > Oy. If not, a jump ahead is executed. If 
it is, then a test is made to determine if the sequence has 
been gone through C times. If so, the program will reini- 
tialize itself by entering the reset loop for a new start. 

This does not put the hybrid system in reset state but only 
restarts the algorithm. If the program did not reinitialize, 
then a is set to a lf to prepare it for another cycle 
through the sequence. 

This brings the execution of the program to the generation 
of the noise vector | k+1 (where it would have been if there 
were less than q consecutive failures or if cr was not 
greater than Oy) . The distribution of the noise generator 
is uniform over the space [-W, W] if no success has been 
achieved since first entering the search loop. Otherwise the 
distribution is gaussian with zero mean. Next, the new 

lr_L. 1 

adjoint initial conditions are obtained by adding | 



to the last set of adjoint initial conditions which produced 
a success, and the program moves to the output section. 

The search mode continues reducing the metric by the 
selection of p(0) until a value is found which gives a 
solution to the simulated control problem. This solution is 
sensed at the point in the program where is compared 

with e. This test results in a branch to the left where 
the first operation is to set up the end -state loop condi- 
tion. Next, the data representing a solution are logged on 
the printer and the new adjoint initial conditions set to 
the value which gave the solution. The program then goes 
to output with the system set to perform the end-state loop. 

C. End -State Loop — The end -state loop is shown in 
figure 8 by the line weight emphasis . In end state the 
integrators still sequence through their initial condition 
or operate modes just as in the search mode. However, 
p(0) remains fixed at that value which produced the solution. 

In this manner it is possible to observe the optimal solutions 

* 

found by the search algorithm on the CRT display described 
earlier . 

AN EXAMPLE ORBIT -TRANSFER PROBLEM 

In this section the usefulness of the random search 
algorithm in solving a moderately difficult problem will be 
discussed. An orbit -transfer problem was selected because 



it is high order, nonlinear, and is the type of problem for 
which optimization is appropriate. That is, since such large 
amounts of fuel are involved in effecting an orbit transfer, 
it would be profitable to minimize the required fuel. In the 
first four subsections, we will formulate the problem, show 
the equations which result from an application of the Maximum 
Principle, and discuss the boundary conditions and the asso- 
ciated metric; in the last four subsections, we will present 
some experimental results on the characteristics of the random 
search algorithm and on the orbit -transfer problem. 

Formulation 

The particular orbit -transfer problem considered was a 
transfer from an earth-moon trajectory to a circular orbit 
around the moon. The physical situation is illustrated in 
figure 9 for "the planar situation. The final desired orbit 
is a circle 190 km above the lunar surface, i.e., a circle 
of radius 1928 km. The midcourse corrections might have 
placed the vehicle on some typical orbits as shown in which 
pericynthian may or may not coincide with the desired final 
orbit. The equations of motion in the vicinity of the moon 
will be governed by two-body equations and the trajectories 
will be hyperbolas as indicated. For purposes of simulation, 
a good coordinate system to describe the motion is the rotating 
coordinate system shown in figure 9 ; the origin is on the 
circular orbit and its velocity is the same as for 
a particle in that circular orbit. This coordinate system 
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is a desirable one to use because it subtracts out large 
constant values from the inertial coordinate system. Assuming 
the vehicle is to be controlled by gimbaling a thrusting 
engine in which the mass flow rate is used to vary the thrust , 
the exact equations of motion in the rotating coordinate 
system shown in figure 9 are: 

m(t)X - 2wm(t)Y = -m(t)c cos a + m(t)(^oj 2 - jt^X 


+ M(t)(w 2 - r 


) (15) 


m(t)Y + 2tom(t)X = -m(t)c cos 3 + m(t)fu) 2 — tHY J 


The gimbaling implies the constraint 


cos 2 a + cos 2 P = 1 


(16) 


In the above equations m(t) is the vehicle mass, u> is the 
angular velocity of the rotating coordinate system, c is 
the exhaust velocity, a and 3 are the angles between the thrust 
vector and the X and ^ axes, respectively, r is the radius 
to the vehicle in the inertial coordinate system, r c is the 
radius of the desired circular orbit, and p is a gravita- 
tional constant. Typical values used were m(0) = 39*096 kg, 
c = 3.1A05 km/sec, r c = 1928.68 km, p = h,890 km 3 /sec 2 , 
m(t) = -31>07 kg/sec, u = 8.259X10 -4 r/sec. 

Now if we let x x = X, x 2 = X, X3 = Y, x 4 = Y, and 
x 5 = m(t), the equations of motion can be put in the state 
vector form 


x = f(x, u) 



where x = (x x , x 2 , x 3 , x 4 , x 5 ) and u = (u x , u 2 , U3) • In 
expanded form we have the set of five nonlinear differential 
equat ions : 


Xl 

x 2 

X3 

x 4 


= x 2 


„ CUU.U3 

2wx 4 + + 


x 5 


(f - ^) xi + - 


is F c 


= X 4 


) (17) 


CU 2 U 3 / o U ' 

-2wx 2 + -3^ + [J* - zk) x 3 


X 5 = -U3 


where u x = cos a, u 2 = cos 3, and U3 = -m(t) . It is worth 
noting with respect to simulation accuracy that by using the 
rotating coordinate system the terms involving r which are 
difficult to simulate accurately provide small corrections 
to the more dominant terms involving only the states in the 
rotating coordinate system. These correction terms cannot 
be ignored, however- A simpler model which ignores these 
terms was found to introduce significant errors in end 
states and fuel required. 

The task to be accomplished is to minimize the fuel to 
go from a given position on the hyperbolic orbit to the 
circular orbit. Thus, the quantity of fuel used is intro- 
duced as the added coordinate: 

x c (t) = J U3(r)dT (18) 

o 

and we can interpret the objective as the minimization of 
the terminal value x 0 (T) . This additional coordinate 
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requires the corresponding differential equation 

X 0 (t) = U3(t) (19) 

to be adjoined to the set (17) • 


Maximum Principle Solution 


Although application of the Maximum Principle to the 
present problem will not be given in detail here, the equa- 
tions important to the hybrid computer implementation will 
be summarized. First are the adjoint equations which are 
fundamental to the development. The approximate adjoint 
equations can be shown to be 



P 2 “ -Pi + 2wp 4 





( 20 ) 


Second are the equations for the optimal control vector which 
have been derived from the Maximum Principle: 
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P2 






( 21 ) 


U3 = M if v - ||P|| - ^ (P 5 + 1) > 0 

0 otherwise ; 

where |jp|| = J p g 2 + p 4 2 and. M is the magnitude of the 
maximum thrust. It is seen that the thrust magnitude is 
either on or off depending on the sign of the switching 
function v, and the thrust angles are continuous functions 
of the adjoint variables. To obtain an explicit solution, 
it is now necessary to solve simultaneously by means of the 
analog computer the equations ( 17 ) for x, equations (20) 
for p, and equations (21) for control u, subject to 
certain boundary conditions yet to be discussed. The auxil- 
iary differential equations for x c (t) and p 0 (t) do not 
couple to the other equations; they can, therefore, be dropped 
in the implementation to be described later. 

Boundary Conditions 

The boundary conditions, as yet unspecified, are the 
heart of the random search method. The initial values for 
the state vector x(0) are fixed at the given values while 
the initial values of the p(0) vector are chosen by the 
algorithm. The desired terminal conditions can be specified 
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in a variety of ways depending on the target set chosen. 

On the one hand, the end points could be treated as fixed at 
some point such as the origin of the coordinate system. This 
is somewhat restrictive. On the other hand, the end points 
could be treated as variable by taking the target set to be 
the desired circular orbit . This approach would yield the 
best point in the target set but would require complicated 
transversality conditions to be satisfied. We will take 
a middle ground in the interests of not unduly complicating 
the problem. We will take the target set to be the desired 
circular orbit, but we will consider the end point fixed at 
whatever point in the target set is reached in the fixed 
time, T. The advantage of this specification is that the 
boundary conditions are somewhat simpler than if transver- 
sality conditions had been included. With these boundary 
conditions the vehicle will reach some point on the desired 
circular orbit but perhaps not the best point. However, as 
we will see in a later example, the fuels for the large 
majority of solutions were so close to the theoretical mini- 
mum fuel based on impulsive orbit transfer that it does not 
appear that satisfying transversality conditions could result 
in a better solution with lower fuel. This conclusion was 
found valid over the range for which the problem was studied. 



The specific values for the terminal boundary conditions can 
be found from conventional theory and summarized as: 


x x (T), x 2 (T), x 3 (T), x 4 (t) 

fixed 

x 5 (T) 

free 

Pi(T), p 2 (T), p 3 (T), p 4 (T) 

free 


( 22 ) 


P S (T) = 0 fixed J 

The Vector Metric 

For purposes of satisfying the boundary conditions, the 
components of the vector metric J need to be specified. 
Since the target set is taken to be the circular orbit, a 
suitable displacement metric is 


J D - | r(t) - r c 


( 23 ) 


where rc is the radius of the desired circular orbit . As 
for the velocity metric, it is clear that for the vehicle to 
stay close to the desired circular orbit after the terminal 
time T, we would like: (a) the radial component of the 

inertial velocity, V R , to be small, and (b) the tangential 
component, Vrp , to be close to the velocity consistent with 
that for the circular orbit, i.e., r c u>. Thus, a reasonable 
metric which measures the errors in these velocity components 
is 


Jy ~ J( r c w “ Vip) 2 ' + 


Vr c 


( 24 ) 


The quantities Vqi and Vg are simple transformations of the 
states in the rotating coordinate system. The only adjoint 
variable which must be fixed at the terminal end is p 5 so 
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that we use simply 

J p = p 5 2 (t) (25) 


The values to which the components of J must be reduced 
are somewhat dependent on the mission and accuracy require- 
ments after the terminal time. For most of the study values 
were chosen to be 


Gpj = 8 km 

e v = 15 m/s 
£p = 0*5 VOlt 


(26) 


This velocity requirement reduces the terminal inertial veloc - 
ity to less than 1# of the initial value. The adjoint variable 
requirement was readily found experimentally. As mentioned 
in the discussion following equation ( 7 ), values of ep 
several times this value were found to be satisfactory. 

Behavior and Characteristics of the Algorithm 

The implementation of the random search algorithm for 
the orbit -transfer problem followed closely the discussion 
in a previous section and no further details will be given 
here. With this implementation the algorithm was studied 
to determine some of its important properties and the effect 
of some of its possible variations. The results of such 


studies will be discussed in this section. 
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In order to study the algorithm, we will confine our 
interest to one particular situation of the orbit -transfer 
problem. This situation was chosen such that the vehicle is 
approaching the moon on the hyperbolic orbit number 2 as shown 
in figure 9* Pericynthian coincides with the desired final 
circular orbit, which is 190 km above the lunar surface. 

Rather arbitrarily the initial starting point was chosen 
as -37°) and the time allowed for the vehicle maneuver was 
taken to be 600 seconds. For comparison, the time to reach 
pericynthian with zero control is 534 seconds. 

There are two displays which well illustrate the manner 
in which the system iterates and searches for a solution. The 
first display illustrated in figure 10 gives information about 
the successive improved iterations . In this figure is shown 
the successive improved values of the initial adjoint vector 
p(0) and the corresponding decrease in the boundary cost 
functions Jp, Jy, and Jp. It is worth noting the creeping 
nature of the P 3 (0) value to a solution value of 24 volts 
even though the initial square distribution was ±15 volts 
and the maximum standard deviation of the gaussian distribu- 
tion was Gy - 10 volts. This situation occurred often. 

The second display, illustrated in figure 11, provides infor- 
mation about every iteration. Here the values of the boundary 
cost function Jp, Jy, and Jp are shown for every iteration 
and the successful iterations. It can be noted that when 
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no improvement is obtained in the J vector, the search 
gradually enlarges until an improvement is found. At this 
point the search is localized. With this type of display 
it is easy to select a rational value for the lower limit 
of the standard deviation, Ci . 

The convergence time to obtain a solution is certainly 
one of the central considerations in the random search method. 
Since this time is a random variable, data were taken to 
give suitable averages and some indication of their accuracy. 
In terms of the number of iterations, N, required for a solu- 
tion, it was found from 70 solutions for the particular 
situation described above that 

N = 7,72k 
% = 5,142 

Thus, the average length of time to obtain a solution is 
about 1.25 minutes. The value of cr w gives some indication 
of the accuracy to be expected for measurements of other 
situations. For example, if 25 trials are made for another 
situation, we would expect the standard deviation of the 
average to be ajj = a N /*/25 w 1,000. This relation is only 
an indication, however, since the standard deviation, a^, 
would not likely remain the same for other situations. 

Next, we will examine some of the possible variations 
in the algorithm and its parameters to indicate their relative 
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effects on the convergence. Although these effects are studied 
for the given problem situation, the trends are generally valid 
for other situations. 

First, the effect of the threshold strategy is shown in 
figure 12. The value of t) = 0 corresponds to removing the 
threshold strategy. Also, a standard deviation is shown on 
the curve. It is worth noting that some improvement can be 
obtained by increasing t] . 

Second, the effect of the one-step strategy, the end- 
search strategy, and an increased initial search size are 
given below with standard deviations as indicated. It is 


Average number 
of iterations 


Algorithm with: one -step strategy 


end-search strategy 

7,724 ± 610 

normal initial search size 


Algorithm without one- step strategy 

12,707 ± 1,150 

Algorithm without end-search strategy 

7,065 ± 1,626 

Algorithm with twice initial search size 

9,971 ± 766 


seen that the one-step strategy is fairly effective, since 
its deletion increases convergence time by 65$. It was always 
beneficial in terms of convergence, costs so little in search 
time, and was, therefore, a lways employed. On the other hand, 
it is seen there that the end-search strategy is slightly 




deleterious for the given situation and increases convergence 
time by 10$. This increase is not too significant, however, 
in view of the likely standard deviation. _ The end-search 
strategy appears to be of greatest value when the algorithm 
has iterated to the neighborhood of the minimum and has dif- 
ficulty in meeting the boundary conditions; such a situation 
might exist perhaps in the neighborhood of a very sharp 
valley. In these cases, the end -search strategy greatly 
reduces the size of the search and enhances the certainty 
of finding the minimum. The effect of initial search size 
on convergence is interesting. It is surprising to note that 
an increase in the initial search size of 100$ in all five 
adjoint variables (32 times larger volume) resulted in an 
increase of only 30$ in convergence time. 

The initial starting value J° is another parameter in 
the algorithm which is somewhat arbitrary. The effect of 
different starting values, in figure 13, shows that convergence 
is fairly flat within a wide range. At the lower values, the 
convergence time increases sharply because the search approaches 
the pure random search. 

Results for Orbit -Transfer Problem 

In this section we will illustrate in some detail the type 
of results obtained for one particular" situation, the same 
situation as in the preceding section. 
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A typical solution which was obtained by the random 
search algorithm is illustrated in figure 14 by the rather <^Fig^. l 4 
complete set of time histories; shown are the system states 
in the rotating and inertial coordinate systems, the control 
vector, and the adjoint variables. Perhaps the most interest- 
ing are the velocities X2 and X4 in the rotating frame which 
can be seen to be reduced from large values of around -900 m/s 
and 400 m/s to small values in order to reduce the inertial 
velocity components, Vr and Vrp, to near their desired values. 

Also the quantity r - r c , the difference between the vehicle 
terminal radius and the desired circular orbit radius can be 
seen to have been reduced to a small value so that the vehicle 
will end up close to the desired orbit. The terminal inertial 
value ^ is slightly positive indicating the desired orbit 
was attained slightly past pericynthian. The optimal control 
vector is also of interest; it is indicated that the thrust 
should be turned on at about 260 seconds before the fixed 
final time and directed as shown. 

Average performance and its variation are perhaps more 
significant than an individual solution. Observation of 70 
solutions revealed that their time histories were only 
slightly different. In terms of the end states and fuel 
used, the following statistics were obtained: 
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Ave fuel 

= 10,230 

±80 kg 

Ave X. 

= 1,935 

±5 km 

Ave ^ 

= 33 

±8 km 

Ave r(T) - r c 

= 6 

±5 km 

Ave Vip - r c w 

= 2 

±7 m/s 

Ave V R 

= 1 

±12 m/ 


It is worth noting here (fromX , and r - r c ) that the 
set of reachable points is centered slightly further out 
than the desired circular orbit and slightly past pericynthian; 
the standard deviation shows this set is confined to a very 
small region. The inertial velocities indicate that they have 
been substantially reduced from the initial total magnitude 
of approximately 2400 m/s. The fuel data are interesting 
because the figures are so close to the idealized minimum- 
fuel impulsive orbit transfer of 10,120 kg. 

It is significant to note that from this same set of 
70 solutions, the solutions for the initial condition adjoint 
vector p(0) were vastly different even though the time 
histories of the system states Were so close. The extremes 
observed as well as the p(0) vector corresponding to fig- 
ure 14 are as follows . 


Adjoint 

variable 

Extremes 

Values for 
figure Ik 

Pi 

-31.98 to +7.15 

-11.36 

P 2 

+1.08 to +14.67 

+8.96 

P 3 

-8.33 to +35-60 

-54 

P 4 

-11.94 to + 5-45 

- 9-25 

P* 

-7.98 to +8.11 

-6.30 
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These results are significant because they illuminate the nature 
of the mapping p(0) J. In geometrical terms, discussed in 
more detail later, these results mean that the boundary cost func- 
tion hypersurfaces have many stalactites protruding below the 
€ level. 

Variations in the time allowed to reach the desired orbit 
were studied with much the same type of results as in the 
preceding two paragraphs. Solutions were obtained for times 
from 550 seconds to 85O seconds (the time to reach pericynthian 
with zero control was 534 seconds). The main difference in 
results obtained was that as the allowable time increases, the 
set of reachable points as given by x and ^ moves further 
around the desired circular orbit and the fuel requirements 
appear to increase. For example, for T = 850 seconds, 

1828 km, « 658 km, and the fuel required increased by 
about 300 kg. When the allowable time is outside the range 
given above, no solutions could be obtained. 


Results for Other Situations 


The algorithm was experimentally tested on a variety of 
situations. They include different initial points along each 
of the three trajectories as indicated by the grid of points 
in figure 9, different allowable times for the vehicle to 
reach the desired point, and variations in some of the 
vehicle design parameters. A sampling of some of these data 
and a few comments are appropriate. 

Figures 15 and l6 show two different solutions to the < 
same problem in which the starting point is further along on 


Figs. 15 
ind 16 
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orbit 2 (-13°) than in the previous two sections- These 
solutions are interesting because of the unusual control 
functions found. One of the control solutions has two thrust- 
ing periods (on-off -on) with the initial thrusting being quite 
long; the other solution has only one thrusting period (on-off) 
and again the initial thrusting is long. In contrast, most 
of the solutions found for other situations have only a single 
thrusting period of the off -on type. 

Next, figures 17 and 18 illustrate solutions which were /Figs. 17 

^\^nd l8 


obtained on orbits 1 and 3* On orbit 1, the solutions showed 
there was always one thrusting period of the off -on type (see 
fig. 17) for all starting points along the orbit and all times 
T for which solutions could be obtained; the fuels used were 
close to the values obtained on orbit 2. On orbit 3> the 
solutions nearly always consisted of the more unusual two 
thrusting periods of the on-off-on type (see fig. 18) . It 
was also noted that solutions were quite difficult to obtain 
on orbit 3 and that the fuels required were approximately 
1500 kg greater than on orbits 1 and 2. No explanation was 
apparent . 

Finally, it was demonstrated that the random search 
approach could be of considerable value in preliminary vehicle 
design. For this purpose the effects of varying some of the 
initial vehicle design parameters were investigated, namely, 
thrust level and mass. It is significant to note that solutions 
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could readily be obtained with engines of 75$ and 50$ of 
the hormal thrust capability used in the previous portions of 
the study. Similarly, there were no difficulties in obtaining 
solutions for different initial vehicle masses. 


Boundary Cost Function Surfaces 

Perhaps the most illuminating results were revealed by an 
experimental study of the mapping: p(0) -» J. Geometrically 

this mapping can be interpreted as three hypersurfaces (rep- 
resenting the components of the metric J), which are functions 
of the five adjoint initial conditions. The character of the 


cross sections through these surfaces at a solution point is 
given in figures 19 and 20 for the two different problem 
situations indicated. These curves were obtained by slowly 


varying one of the adjoint variables from -100 volts to +100 volts 
while all the others were fixed at their respective solution 
values. In this manner cross sections through the hypersurfaces 
passing through a solution point were plotted. Note the dip 


in all three surfaces at the optimum value of the adjoint 


variable. The interest in these curves is, of course, in 
the rather irregular, mult i -valley ed nature of the surfaces 
as well as the rather narrow valleys surrounding the optimum 
point . The sharpness or narrowness of this valley gives an 
indication of the difficulty in finding the solution. Further, 
these surfaces clearly indicate the difficulties which gradient 
methods would have because of the likely possibility of "hanging 
up" in the wrong valley. The two-dimensional cross sections 



shown give only a hint of the difficulties to be encountered 
in the actual six -dimensional space. 

Another more pictorial representation (for the same 
problem situation as in figure 19) is given by the three- 
dimensional views in figure 21. Here are shown the boundary 
cost function surfaces in the pi, P2 plane and in the p 1? 
p 4 plane while all other adjoint variables are fixed at their 
solution values. These results illustrate the character of 
the surfaces away from the optimum point in two directions 
instead of only one direction. 

Further clues as to the nature of the boundary cost 
function surfaces are revealed by previously presented data 
in which it was shown that many distinctly different p(0) 
vectors result in solutions. The significance of this type 
data is that it shows that the surfaces have many stalactites 
which protrude below the e level. 

CONCLUDING REMARKS 

It is worth emphasizing that the objective in the' present 
study has been to devise and demonstrate a feasible method of 
implementing the Maximum Principle based on an adaptive random 
search algorithm. It has been demonstrated that the computa- 
tional requirements are most effectively implemented on a hybrid 
computer system. The method has been shown to be a practical 
approach to generating explicit optimal solutions for a moderately 
complex example problem under a wide range of situations. The 
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versatility and ease in obtaining solutions makes it a valuable 
approach for preliminary system design because it allows one 
to study the tradeoffs among various initial design choices in 
terms of performance. More generally, the method could be 
applied to any two -point boundary value problem or to parameter 
optimization of systems with given configurations. 

More effective utilization of the method either in an 
on-line or off-line implementation could be realized by an 
increase in the iteration rate. Such increases will in 
general have to be the result of advances in hybrid systems . 
However, a moderate increase in rate could be realized with 
the present system by careful redesign of the program. 

Undoubtedly the most important limitation of the method 
is due to the minimum resolution of the analog computer. 

Thus certain problems with excessive dynamic range might be 
excluded. 

The approach has a number of advantages over alternate 
approaches. It is conceptually simple, straightforward, and 
easily implemented. It can be applied uniformly within its 
limitations. No restrictions are placed on the boundary 
surface, and an analytical expression for the surface is not 
required. The algorithm's local and global search properties 
provide it with the ability to jump from one local valley to 
another so that "hanging up" in local valleys is avoided. 
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FIGURE LEGENDS 

Figure 1.- Hybrid system block diagram for random search method. 
Figure 2.- Typical boundary cost function surface. 

Figure 3*- Hybrid system hardware. 

Figure 4.- Sequencing of events during one iteration. 

Figure 5* - Algorithm flow graph. 

Figure 6.- Reset loop. 

Figure 7*- Search loop. 

Figure 8.- End-state loop. 

Figure 9-- Geometry of orbit -transfer problem. 

Figure 10.- Behavior of successful iterations during search. 
Figure 11.- Behavior of every iteration during search. 

Figure 12.- Effect on convergence of threshold strategy. 

Figure 13 • - Effect on convergence of varying initial value J°. 
Figure 14.- Time history solutions; orbit 2, x(0) at -37°* 

T = 600 s ec . 

Figure 15*- Time history solutions 1; orbit 2 y x(0) at -13°, 

T - 400 sec . 

Figure l6.- Time history solutions 2; orbit 2, x(0) at -13°, 

T - 400 sec . 

Figure 17.- Time history solutions; orbit 1, x(0) at -6o°, 

T =• 971 sec. 

Figure l8.- Time history solutions; orbit 3> x(0) at -60°, 


T = 1075 sec. 
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Figure 19-- Two-dimensional cross sections of boundary hyper- 
surface; orbit 2, x(0) at - 31 ° } T = 600 sec. 

Figure 20.- Two-dimensional cross sections of boundary hyper - 
surface; orbit 3 , x(0) at - 60 ° , T = 1075 sec. 
Figure 21.- Three-dimensional cross sections of boundary hyper 
surface; orbit 2 , x(0) at - 31 ° > T = 600 sec. 

(a) px , p 2 plane. 

Figure 21.- Concluded. 

(b) Pi, P4 plane. 




Figure 1.- Hybrid system block diagram for random search method. 









£2 

i\ 

2 § 

U 

U 


< “ 

l| 

il 


r NoiioNnj isoo AavaNnoa 


Figure 2.- Typical boundary cost function surface. 



Figure 3 • - Hybrid system hardware 
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Figure 5 -- Algorithm flow graph. 
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Figure 6.- Reset loop. 
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Figure 7*" Search loop. 
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Figure 8.- End-state loop. 
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Figure IT-- Time history solutions; orbit 1 , x(0) at -60 ° , T = 971 sec - 
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Figure l6.- Time history solutions 2; orbit 2, x(0) at -13°> T 1+00 sec 
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(a) pi, p 2 plane. 

Figure 21.- Three-dimensional cross sections of boundary hypersurface; 

orbit 2, x(0) at -37°, T = 600 sec . 
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(b) pi, P 4 plane . 
Figure 21.- Concluded. 
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